// Copyright 2020 Poofee (https://github.com/Poofee)
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
// ------------------------------------------------------------------------
/*****************************************************************************
 *                                                                           *
 *  File:    relay2d.cpp                                                     *
 *                                                                           *
 *****************************************************************************
 *                                                                           *
 *  Authors: Poofee                                                          *
 *  Email:   poofee@qq.com                                                   *
 *  Address: https://github.com/Poofee                                       *
 *  Date:    2020年07月15日                                                   *
 *                                                                           *
 *****************************************************************************/
#include "relay2d.h"
#include "SuperLU_MT.h"

Relay2D::Relay2D() {

}



/*!
 * @brief      二维重分网。
 * @author     Poofee
 * @date       2020年9月18日
 * @param[out] 
 * @param[in]  int timestep  
 * @return     void  
**/
void Relay2D::remesh(int timestep) {
	/** 单位换算，不然不对 **/
	Relay::remesh(0, Ddisplacements[timestep] / unitRatio);
}


/*!
 * @brief      计算三角形单元的基本项
 * @author     Poofee
 * @date       2020年9月18日
 * @param[out] 
 * @param[in]  int index  
 * @return     void  
**/
void Relay2D::makeTriangle(int index) {
    int k,m,n;
    double p0,p1,p2,q0,q1,q2,area;
    k = pmeshele[index].n[0];
    m = pmeshele[index].n[1];
    n = pmeshele[index].n[2];

    p0 = pmeshnode[m].y - pmeshnode[n].y;
    pmeshele[index].P[0] = p0;
    p1 = pmeshnode[n].y - pmeshnode[k].y;
    pmeshele[index].P[1] = p1;
    p2 = pmeshnode[k].y - pmeshnode[m].y;
    pmeshele[index].P[2] = p2;

    q0 = pmeshnode[n].x - pmeshnode[m].x;
    pmeshele[index].Q[0] = q0;
    q1 = pmeshnode[k].x - pmeshnode[n].x;
    pmeshele[index].Q[1] = q1;
    q2 = pmeshnode[m].x - pmeshnode[k].x;
    pmeshele[index].Q[2] = q2;

    area = 0.5*std::abs(p1 * q2 - q1 * p2);
    pmeshele[index].AREA = area;
    pmeshele[index].rc = (pmeshnode[k].x +
                          pmeshnode[m].x +
                          pmeshnode[n].x) / 3;
    pmeshele[index].zc = (pmeshnode[k].y +
                          pmeshnode[m].y +
                          pmeshnode[n].y) / 3;

    int flag = 0;
    for (int f = 0; f < 3; f++) {
        if (pmeshnode[pmeshele[index].n[f]].x < 1e-7) {
            flag++;
        }
    }
    /** 计算三角形重心半径 **/
    if (flag == 2) {
        pmeshele[index].ydot = pmeshele[index].rc;
    } else {
        pmeshele[index].ydot  = 1 / (pmeshnode[k].x + pmeshnode[m].x);
        pmeshele[index].ydot += 1 / (pmeshnode[k].x + pmeshnode[n].x);
        pmeshele[index].ydot += 1 / (pmeshnode[m].x + pmeshnode[n].x);
        pmeshele[index].ydot = 1.5 / pmeshele[index].ydot;
    }

//    double Y11,Y12,Y13,Y22,Y23,Y33;
    pmeshele[index].Y11 = q0 * q0 + p0 * p0;
    pmeshele[index].Y12 = q0 * q1 + p0 * p1;
    pmeshele[index].Y13 = q0 * q2 + p0 * p2;
    pmeshele[index].Y22 = q1 * q1 + p1 * p1;
    pmeshele[index].Y23 = q1 * q2 + p1 * p2;
    pmeshele[index].Y33 = q2 * q2 + p2 * p2;

    pmeshele[index].Y11 /= 4. * area;
    pmeshele[index].Y12 /= 4. * area;
    pmeshele[index].Y13 /= 4. * area;
    pmeshele[index].Y22 /= 4. * area;
    pmeshele[index].Y23 /= 4. * area;
    pmeshele[index].Y33 /= 4. * area;
}


/*!
 * @brief      
 * @author     Poofee
 * @date       2020年9月18日
 * @param[out] 
 * @return     void  
**/
void Relay2D::setBoundary() {
	/** 初始化，提出读取的gmsh多余的几何节点以及边界点 **/
	int node_bdr = boundaryPoints.size();
	int node_all = allPoints.size();

	for (int i2 = 0; i2 < num_pts; ++i2) {
		pmeshnode[i2].bdr = 3;
	}
	for (int i2 = 0; i2 < node_all; ++i2) {
		pmeshnode[allPoints[i2]].bdr = 0;
	}
	for (int i2 = 0; i2 < node_bdr; ++i2) {
		pmeshnode[boundaryPoints[i2]].bdr = 3;
	}
	/** 将边界点排序到末尾 **/
	node_bdr = 0;
	for (int ibdr = 0; ibdr < num_pts; ++ibdr) {
		if (pmeshnode[ibdr].bdr == 3) {
			++node_bdr;
			node_reorder(num_pts - node_bdr) = ibdr;
			node_pos(ibdr) = num_pts - node_bdr;
			pmeshnode[ibdr].A = 0;
			A(ibdr) = 0;
		} else {
			node_reorder(ibdr - node_bdr) = ibdr;
			node_pos(ibdr) = ibdr - node_bdr;
		}
	}
}


/*!
 * @brief      
 * @author     Poofee
 * @date       2020年9月18日
 * @param[out] 
 * @return     void  
**/
void Relay2D::allocateMemory() {
	int size = 4;
	locs.resize(2, size * size * num_triangle);
	vals.resize(1, size * size * num_triangle);
	bbJz.resize(num_pts); bbJz.zeros();
	node_reorder.resize(num_pts);
	node_pos.resize(num_pts);
	bn.resize(num_pts);
	A.resize(num_pts); A.zeros();
    //A_old.resize(num_pts);/** 上一步的大小肯定和这一步不一样 **/
}


/*!
 * @brief      二维轴对称静态
 * @author     Poofee
 * @date       2020年9月18日
 * @param[out] 
 * @return     void  
**/
void Relay2D::MagnetoStatic2DAxisNR() {
    /** 计时 **/
    clock_t time[10];
    int tt = 0;
    time[++tt] = SuperLU_timer_();

    double ce[3][3] = { {0} };/** 雅可比矩阵 **/
    double Je[3][3] = { {0} };/** 用于牛顿迭代 **/
    double De[3]    = {  0  };/** 电流密度矩阵 **/

    double Y11, Y12, Y13, Y22, Y23, Y33;
    /** 时间步循环 **/
    current_step = 1;
    double error;
    int timesize = timesteps.size();
    for (int i = 1; i < timesize; ++i) {
        printCurrentTime();
        printf("Time step %d/%d.\n", i, timesize);
        /** remesh，要使用增量位移，向下为负，向上为正 **/
        remesh(i);
        /** 读取第i步的分网 **/
        loadMesh();
        unitConvert(unitRatio);
        /** 查找边界，默认一次边界 **/
        findBoundaryEdges2D(-1);
        findBoundaryPoints2D(-1);

        A_last = A;
        A_old = A;

        allocateMemory();
        setBoundary();
        findAir();
        initMu();
        /** 自由度：未知节点+ **/
        int num_dof = allPoints.size() - boundaryPoints.size();
        printf("DoF %d\n", num_dof);
        double* unknown_b = (double*)calloc(num_dof, sizeof(double));
        double miuold = 1;

        for (int i_tri = 0; i_tri < num_triangle; ++i_tri) {
            int i1 = num_ele - num_triangle + i_tri;
            CElement* ele = pmeshele + i1;
            CMaterial* mat = getMaterial(ele->geometry_tag);
            double dtmp = mat->tau * ele->AREA / 3;
            for (int j = 0; j < 3; ++j) {
                bbJz(ele->n[j]) += dtmp;
            }
        }

        /** 非线性迭代求解 **/
        for (int non_iter = 0; non_iter < MAX_NONLINEARSTEPS; ++non_iter) {
            printCurrentTime();
            printf("Newton's step %d(%d/%d)\n", i, non_iter, MAX_NONLINEARSTEPS);
            /** 有限元装配 **/
            int pos = 0;
            bn.zeros();
            for (int i_tri = 0; i_tri < num_triangle; ++i_tri) {
                int k, m, n;
                int i1 = num_ele - num_triangle + i_tri;
                CElement* ele = pmeshele + i1;
                k = ele->n[0]; m = ele->n[1]; n = ele->n[2];

                /** 计算除了磁导率部分的系数 **/
                makeTriangle(i1);

                /** 计算雅可比矩阵 **/
                double ydot = ele->ydot;
                double miut = ele->miut;
                miuold = ele->miu;
                /** 相当于C矩阵 **/
                Y11 = ele->Y11; Y22 = ele->Y22; Y33 = ele->Y33;
                Y12 = ele->Y12; Y13 = ele->Y13; Y23 = ele->Y23;

                /** 应当是当前迭代的mu **/
                ce[0][0] = Y11;      ce[1][1] = Y22;      ce[2][2] = Y33;
                ce[0][1] = Y12;      ce[0][2] = Y13;      ce[1][2] = Y23;
                ce[1][0] = ce[0][1]; ce[2][0] = ce[0][2]; ce[2][1] = ce[1][2];

                double v[3];
                v[0] = Y11 * A(k) + Y12 * A(m) + Y13 * A(n);
                v[1] = Y12 * A(k) + Y22 * A(m) + Y23 * A(n);
                v[2] = Y13 * A(k) + Y23 * A(m) + Y33 * A(n);

                CMaterial* mat = getMaterial(ele->geometry_tag);

                double tmp = 0;
                if (non_iter > 0) {
                    if (mat->BHpoints == 0) {
                        tmp = 0;
                    } else {
                        tmp = mat->getdvdB(ele->B);
                        if (ele->B > 1e-9) {
                            tmp /= ele->B * ele->AREA;  // B==0?
                            tmp /= ydot * ydot * ydot;
                        }
                    }
                }

                Je[0][0] = v[0] * v[0] * tmp; Je[1][1] = v[1] * v[1] * tmp; Je[2][2] = v[2] * v[2] * tmp;
                Je[0][1] = v[0] * v[1] * tmp; Je[0][2] = v[0] * v[2] * tmp; Je[1][2] = v[1] * v[2] * tmp;
                Je[1][0] = Je[0][1];          Je[2][0] = Je[0][2];          Je[2][1] = Je[1][2];

                /** 计算右侧向量 **/
                double dtmp = mat->tau * ele->AREA / 3;
                De[0] = dtmp; De[1] = dtmp; De[2] = dtmp;

                for (int kk = 0; kk < 3; ++kk) {
                    for (int mm = 0; mm < 3; ++mm) {
                        ce[kk][mm] /= miut * ydot;
                        Je[kk][mm] += ce[kk][mm];
                    }
                }

                //double jr = ele->AREA*mat->Jr / 3;
                for (int j = 0; j < 3; ++j) {
                    //bbJz(ele->n[j]) += jr;
                    /** 计算永磁部分 **/
                    bn(ele->n[j]) -= mat->H_c / 2. * ele->Q[j];

                    bn(ele->n[j]) -= ce[j][0] * A(k) + ce[j][1] * A(m) + ce[j][2] * A(n);
                }

                int n_row, n_col;
                for (int row = 0; row < 3; ++row) {
                    n_row = node_pos(ele->n[row]);
                    for (int col = 0; col < 3; ++col) {
                        /** 判断节点是否在未知节点内 **/
                        /** 得到排序之后的编号 **/
                        n_col = node_pos(ele->n[col]);
                        if (n_row < num_dof && n_col < num_dof) {
                            insertSparseMatrixData(pos++, n_row, n_col, Je[row][col]);
                        }
                    }
                }
            }/** 单元循环结束 **/
            if (non_iter == 0) {
                locs.reshape(2, pos);
                vals.reshape(1, pos);
            }
            bn += bbJz;
            /** 调用线性求解器求解 **/
            /** 使用构造函数来生成稀疏矩阵 **/
            sp_mat X(true, locs, vals, num_dof, num_dof, true, true);

            for (int idof = 0; idof < num_dof; ++idof) {
                unknown_b[idof] = bn(node_reorder(idof));
            }
            //---------------------superLU_MT---------------------------------------
            CSuperLU_MT superlumt(num_dof, X, unknown_b, ParallelThreads);
            if (superlumt.solve() == 1) {
                printf("Error: superlumt.slove. Info:%d\n", superlumt.info);
                return;
            } else {
                A_old = A;
                double* sol = superlumt.getResult();

                for (int idof = 0; idof < num_dof; ++idof) {
                    pmeshnode[node_reorder(idof)].A = sol[idof];// / pmeshnode[i].x;//the A is r*A_real
                    A(node_reorder(idof)) = sol[idof];
                }
                //printf("norm((A), 2):%lf\n", norm((A), 2));
                error = norm((A), 2) / norm(A_old, 2);
                A += A_old;
            }
            /** 更新磁场结果 **/
            double bx, by;
            for (int i_tri = 0; i_tri < num_triangle; ++i_tri) {
                bx = 0; by = 0;
                CElement* ele = pmeshele + num_ele - num_triangle + i_tri;
                for (int j = 0; j < 3; ++j) {
                    bx += ele->Q[j] * A(ele->n[j]);
                    by -= ele->P[j] * A(ele->n[j]);
                }
                CMaterial* mat = getMaterial(ele->geometry_tag);
                ele->B = sqrt(bx * bx + by * by) / 2. / ele->AREA / ele->ydot;
                ele->Bx = bx / 2. / ele->AREA / ele->ydot;
                ele->By = by / 2. / ele->AREA / ele->ydot;
                ele->miut = mat->getMiu(ele->B);
            }

            printf("error:%lf\n", error);
            if (error < Precision) {
                break;
            }
            bn.zeros();
            pos = 0;
        }/** 牛顿迭代结束 **/

        /** 计算电磁力 **/
        calcMagForce(-1, FORCEMETHOD::VIRTUALWORK);

        for (int i_tri = 0; i_tri < num_triangle; ++i_tri) {
            int i1 = num_ele - num_triangle + i_tri;
            pmeshele[i1].miu = pmeshele[i1].miut;
        }

        char fn[256];
        sprintf(fn, "%s%d%sAnoair.vtk", fileName, i, prefixName);
        writeVtkFileNoAir(fn);
        /** 更新结果 **/
        stepIncrement();
        Ddisplacements[current_step] = 0;
        displacements[current_step]  = 0;
        velocities[current_step]     = 0;
        accelerations[current_step]  = 0;
        PhiCoil[current_step]        = 0;
        ICoil[current_step]          = 0;
        magforces[current_step]      = 0;
        UCoil[current_step]          = 0;

        /** 回收空间 **/
        if (unknown_b) delete unknown_b;
    }
}


/*!
* @brief      使用牛顿迭代法计算瞬态感应现象，只有简单的线圈。变量为i。
* @author     Poofee
* @date       2020年9月2日
* @param[out] 
* @return     void  
**/
void Relay2D::MagnetoDynamics2DAxisNR() {
	/** 计时 **/
	clock_t time[10];
	int tt = 0;
	time[++tt] = SuperLU_timer_();
	/** 电路参数 **/

	double ce[3][3] = { {0} };/** 雅可比矩阵 **/
	double Je[3][3] = { {0} };/** 用于牛顿迭代 **/
	double Te[3][3] = { {0} };/** 涡流矩阵 **/
	double De[3] = { 0 };/** 电流密度矩阵 **/

	double Y11, Y12, Y13, Y22, Y23, Y33;
	/** 变量 **/
	double beta = 1;
	/** 时间步循环 **/
	current_step = 1;
	double error;
	double totalTime = 0;
	int timesize = timesteps.size();
	for (int i = 1; i < timesize; ++i) {
        printCurrentTime();
		double dt = timesteps[i];
		totalTime += dt;
        printf("Time step %d/%d, current time is %lf s.\n", i, timesize, totalTime);
		/** remesh，要使用增量位移，向下为负，向上为正 **/
		remesh(i);
		/** 读取第i步的分网 **/
		loadMesh();
		unitConvert(unitRatio);
		/** 查找边界，默认一次边界 **/
		findBoundaryEdges2D(-1);
		findBoundaryPoints2D(-1);

		//testBoundaryResults((char*)"../geo/JRS1250/bnd.vtk");

		A_last = A;
		A_old = A;
		if (i == 0) {
			A.resize(num_pts);
			A.zeros();
			A_old = A;
			A_last = A;
		} else {

		}
		allocateMemory();
		setBoundary();
		findAir();
		initMu();
		/** 自由度：未知节点+ **/
		int num_dof = allPoints.size() - boundaryPoints.size();
		printf("Dof %d\n", num_dof);
		double icoil = UCoil[i - 1];
		double* unknown_b = (double*)calloc(num_dof, sizeof(double));
		double miuold = 1;

		for (int i_tri = 0; i_tri < num_triangle; ++i_tri) {
			int i1 = num_ele - num_triangle + i_tri;
			CElement* ele = pmeshele + i1;
			CMaterial* mat = getMaterial(ele->geometry_tag);
			double dtmp = mat->tau * ele->AREA / 3;
			for (int j = 0; j < 3; ++j) {
				bbJz(ele->n[j]) += dtmp;
			}
		}

		/** 非线性迭代求解 **/
		for (int non_iter = 0; non_iter < MAX_NONLINEARSTEPS; ++non_iter) {
            printCurrentTime();
			printf("Newton's step %d(%d/%d)\n", i, non_iter, MAX_NONLINEARSTEPS);
			/** 有限元装配 **/
			int pos = 0;
			bn.zeros();
			//bbJz.zeros();
			//unknown_b[num_dof] = 0;
			for (int i_tri = 0; i_tri < num_triangle; ++i_tri) {
				int k, m, n;
				int i1 = num_ele - num_triangle + i_tri;
				CElement* ele = pmeshele + i1;
				k = ele->n[0]; m = ele->n[1]; n = ele->n[2];

				/** 计算除了磁导率部分的系数 **/
				makeTriangle(i1);

				/** 计算雅可比矩阵 **/
				double ydot = ele->ydot;
				double miut = ele->miut;
				miuold = ele->miu;
				/** 相当于C矩阵 **/
				Y11 = ele->Y11; Y22 = ele->Y22; Y33 = ele->Y33;
				Y12 = ele->Y12; Y13 = ele->Y13; Y23 = ele->Y23;

				/** 应当是当前迭代的mu **/
				ce[0][0] = Y11;      ce[1][1] = Y22;      ce[2][2] = Y33;
				ce[0][1] = Y12;      ce[0][2] = Y13;      ce[1][2] = Y23;
				ce[1][0] = ce[0][1]; ce[2][0] = ce[0][2]; ce[2][1] = ce[1][2];

				double v[3];
				v[0] = Y11 * A(k) + Y12 * A(m) + Y13 * A(n);
				v[1] = Y12 * A(k) + Y22 * A(m) + Y23 * A(n);
				v[2] = Y13 * A(k) + Y23 * A(m) + Y33 * A(n);

				CMaterial* mat = getMaterial(ele->geometry_tag);

				double tmp = 0;
				if (non_iter > 0) {
					if (mat->BHpoints == 0) {
						tmp = 0;
					} else {
						tmp = mat->getdvdB(ele->B);
						if (ele->B > 1e-9) {
							tmp /= ele->B * ele->AREA;  // B==0?
							tmp /= ydot * ydot * ydot;
						}
					}
				}

				Je[0][0] = v[0] * v[0] * tmp; Je[1][1] = v[1] * v[1] * tmp; Je[2][2] = v[2] * v[2] * tmp;
				Je[0][1] = v[0] * v[1] * tmp; Je[0][2] = v[0] * v[2] * tmp; Je[1][2] = v[1] * v[2] * tmp;
				Je[1][0] = Je[0][1];          Je[2][0] = Je[0][2];          Je[2][1] = Je[1][2];

				/** 计算右侧向量 **/
				double dtmp = mat->tau * ele->AREA / 3;
				De[0] = dtmp; De[1] = dtmp; De[2] = dtmp;
				double dtmp1 = 0;// 2 * PI * Yc * dtmp;
				/** 计算涡流矩阵 **/
				double etmp = 0;// mat->sigma / ydot * ele->AREA / 12;
				Te[0][0] = etmp * 2 + dtmp1 * bbJz(k); Te[1][1] = etmp * 2 + dtmp1 * bbJz(m); Te[2][2] = etmp * 2 + dtmp1 * bbJz(n);
				Te[0][1] = etmp * 1 + dtmp1 * bbJz(m); Te[0][2] = etmp * 1 + dtmp1 * bbJz(n); Te[1][2] = etmp * 1 + dtmp1 * bbJz(n);
				Te[1][0] = etmp * 1 + dtmp1 * bbJz(k); Te[2][0] = etmp * 1 + dtmp1 * bbJz(k); Te[2][1] = etmp * 1 + dtmp1 * bbJz(m);

				for (int kk = 0; kk < 3; ++kk) {
					for (int mm = 0; mm < 3; ++mm) {
						//Te[kk][mm] /= dt;
						ce[kk][mm] /= miut * ydot;
						Je[kk][mm] += ce[kk][mm];
						//Je[kk][mm] += Te[kk][mm];

						//ce[kk][mm] += Te[kk][mm];
					}
				}

				//double jr = ele->AREA*mat->Jr / 3;
				for (int j = 0; j < 3; ++j) {
					//bbJz(ele->n[j]) += jr;
					/** 计算永磁部分 **/
					bn(ele->n[j]) -= mat->H_c / 2. * ele->Q[j];

					//if (i_tri < air_position) {
					bn(ele->n[j]) -= ce[j][0] * A(k) + ce[j][1] * A(m) + ce[j][2] * A(n);// -Yc * De[j] * icoil;
				//}
				//if (i_tri < air_position && i > 1) {
				//	//if(mat->sigma > 0)
				//	bn(ele->n[j]) += Te[j][0] * A_last(k) + Te[j][1] * A_last(m) + Te[j][2] * A_last(n);
				//}
				}

				int n_row, n_col;
				for (int row = 0; row < 3; ++row) {
					n_row = node_pos(ele->n[row]);
					for (int col = 0; col < 3; ++col) {
						/** 判断节点是否在未知节点内 **/
						/** 得到排序之后的编号 **/
						n_col = node_pos(ele->n[col]);
						if (n_row < num_dof && n_col < num_dof) {
                            insertSparseMatrixData(pos++, n_row, n_col, Je[row][col]);
						}
					}

					//if (n_row < num_dof) {
					//	locs(0, pos) = n_row;
					//	locs(1, pos) = num_dof;
					//	vals(0, pos) = -Yc * De[row];
					//	pos++;
					//	locs(0, pos) = num_dof;
					//	locs(1, pos) = n_row;
					//	vals(0, pos) = -Yc * De[row];
					//	pos++;
					//}
				}
				//if (i_tri < air_position) {
				//	unknown_b[num_dof] -= -Yc * De[0] * A(k) - Yc * De[1] * A(m) - Yc * De[2] * A(n);
				//}
				//if (i_tri < air_position && i > 1) {
				////if (mat->sigma > 0)
				//	unknown_b[num_dof] += -Yc * De[0] * A_last(k) - Yc * De[1] * A_last(m) - Yc * De[2] * A_last(n);
				//}
			}/** 单元循环结束 **/
			/** 右下角的点 **/
			//locs(0, pos) = num_dof;
			//locs(1, pos) = num_dof;
			//vals(0, pos) = c;
			//pos++;
			if (non_iter == 0) {
				locs.reshape(2, pos);
				vals.reshape(1, pos);
			}
			//bn.zeros();
			//bn += bbJz;
			/** 调用线性求解器求解 **/
			/** 使用构造函数来生成稀疏矩阵 **/
			sp_mat X(true, locs, vals, num_dof, num_dof, true, true);

			for (int idof = 0; idof < num_dof; ++idof) {
				unknown_b[idof] = bn(node_reorder(idof));
			}
			/** 对应电压的部分 **/
			//unknown_b[num_dof] += -c * icoil - (C1 + C2) / 2 / PI * UCoil[i - 1];
			//char fn[256];
            //sprintf(fn, "eddy%s%dunknown_b%d.txt", fileName, i, non_iter);
			//FILE *fp1 = fopen(fn, "w");
			//for (int it = 0; it <= num_dof; ++it) {
			//	fprintf(fp1, "%lf\n", unknown_b[it]);
			//}
			//fclose(fp1);
			//printf("unknown_b[num_dof]:%lf\t", unknown_b[num_dof]);
			//---------------------superLU_MT---------------------------------------
            CSuperLU_MT superlumt(num_dof, X, unknown_b, ParallelThreads);
			if (superlumt.solve() == 1) {
				printf("Error: superlumt.slove. Info:%d\n", superlumt.info);
				return;
			} else {
				A_old = A;
				double* sol = superlumt.getResult();

				for (int idof = 0; idof < num_dof; ++idof) {
					pmeshnode[node_reorder(idof)].A = sol[idof];// / pmeshnode[i].x;//the A is r*A_real
					A(node_reorder(idof)) = sol[idof];
				}
				//printf("norm((A), 2):%lf\n", norm((A), 2));
				error = norm((A), 2) / norm(A_old, 2);
				A += A_old;
				//icoil += sol[num_dof];
				//printf("sol[num_dof]:%lf\t", sol[num_dof]);
				/*char fn[256];
				sprintf(fn, "%dsol%d.txt", i, non_iter);
				fp1 = fopen(fn, "w");
				for (int it = 0; it <= num_dof; ++it) {
					fprintf(fp1, "%lf\n", sol[it]);
				}
				fclose(fp1);*/
			}
			/** 更新磁场结果 **/
			double bx, by;
			for (int i_tri = 0; i_tri < num_triangle; ++i_tri) {
				bx = 0; by = 0;
				int i1 = num_ele - num_triangle + i_tri;
				CElement* ele = pmeshele + i1;
				for (int j = 0; j < 3; ++j) {
					bx += ele->Q[j] * A(ele->n[j]);
					by -= ele->P[j] * A(ele->n[j]);
				}
				CMaterial* mat = getMaterial(ele->geometry_tag);
				ele->B = sqrt(bx * bx + by * by) / 2. / ele->AREA / ele->ydot;
				ele->Bx = bx / 2. / ele->AREA / ele->ydot;
				ele->By = by / 2. / ele->AREA / ele->ydot;
				ele->miut = mat->getMiu(ele->B);

				//printf("%lf\t", ele->B);
			}
			//            

			printf("error:%lf\ticoil:%lf\n", error, icoil);
			if (error < Precision) {
				break;
			}
			bn.zeros();
			bbJz.zeros();
			pos = 0;
		}/** 牛顿迭代结束 **/

		/** 计算电磁力 **/
        calcMagForce(-1, FORCEMETHOD::VIRTUALWORK);

		for (int i_tri = 0; i_tri < num_triangle; ++i_tri) {
			int i1 = num_ele - num_triangle + i_tri;
			pmeshele[i1].miu = pmeshele[i1].miut;
		}


		char fn[256];
        sprintf(fn, "%s%d%sAnoair.vtk", fileName, i, prefixName);
        writeVtkFileNoAir(fn);
		/** 更新结果 **/
		stepIncrement();
		Ddisplacements[current_step] = 0;
		displacements[current_step] = 0;
		velocities[current_step] = 0;
		accelerations[current_step] = 0;
		PhiCoil[current_step] = 0;
		ICoil[current_step] = 0;
		magforces[current_step] = 0;
		UCoil[current_step] = icoil;

		/** 回收空间 **/
		if (unknown_b) delete unknown_b;
	}
}


/*!
 * @brief      
 * @author     Poofee
 * @date       2020年9月18日
 * @param[out] 
 * @return     void  
**/
void Relay2D::findAir() {
	/** 计算固定区域的单元节点数目，起始编号为0 **/
	int i1;
	for (int i_tri = 0; i_tri < num_triangle; ++i_tri) {
		i1 = num_ele - num_triangle + i_tri;
		if (pmeshele[i1].geometry_tag == tag_air) {
			air_position = i_tri;
			break;
		}
	}
}


/*!
 * @brief      
 * @author     Poofee
 * @date       2020年9月18日
 * @param[out] 
 * @return     void  
**/
void Relay2D::initMu() {
	/** 单元初始化，如果是非线性区域，可以考虑优化先初始值设置 **/
	for (int i_tri = 0; i_tri < num_triangle; ++i_tri) {
		int i1 = num_ele - num_triangle + i_tri;
		pmeshele[i1].B = 0;
		/** 对于涡流区域，没有形变，直接将miu带过来 **/
		if (i_tri < air_position && pmeshelelast) {
			pmeshele[i1].miu = pmeshelelast[i1].miut;
			pmeshele[i1].miut = pmeshelelast[i1].miut;
		} else {
			pmeshele[i1].miu = miu0;
			pmeshele[i1].miut = miu0;
		}
    }
}


/*!
 * @brief      写二维的vtk文件。
 * @author     Poofee
 * @date       2020年9月18日
 * @param[out] 
 * @param[in]  char * fname  
 * @return     void  
**/
void Relay2D::writeVtkFile(char* fname) {
	int num_cell = num_triangle;
	int ele_size = 3;
	int ele_type = 5;
	Relay::writeVtkFile(fname, num_cell, ele_type, ele_size);

	FILE* fp;
	fp = fopen(fname, "a+");
	/** 数据 **/
	fprintf(fp, "\nPOINT_DATA %d\n", num_pts);
	fprintf(fp, "SCALARS %s double %d\n", "Az/(Wb/m)", 1);
	fprintf(fp, "LOOKUP_TABLE %s\n", "Atable");
	for (int i = 0; i < num_pts; ++i) {
		if (pmeshnode[i].x > 0) {
			fprintf(fp, "%lf\n", A(i) / pmeshnode[i].x);
		} else {
			fprintf(fp, "%lf\n", A(i));
		}
	}

	fprintf(fp, "\nCELL_DATA %d\n", num_cell);
	fprintf(fp, "SCALARS %s double %d\n", "B/T", 1);
	fprintf(fp, "LOOKUP_TABLE %s\n", "Btable");
	for (int i_tri = 0; i_tri < num_cell; ++i_tri) {
		int i1 = num_ele - num_cell + i_tri;
		fprintf(fp, "%10.8e\n", pmeshele[i1].B);
	}

	fprintf(fp, "\nVECTORS %s double\n","Bvector");
	for (int i_tri = 0; i_tri < num_cell; ++i_tri) {
		int i1 = num_ele - num_cell + i_tri;
        fprintf(fp, "%10.8e %10.8e %10.8e\n", pmeshele[i1].Bx, pmeshele[i1].By, 0.0);
	}

    fclose(fp);
}


/*!
* @brief      目前的程序判断需要把空气保存在后尾，才能避免。
* @author     Poofee
* @date       2020年9月9日
* @param[out] 
* @param[in]  char * fname  
* @return     void  
**/
void Relay2D::writeVtkFileNoAir(char * fname) {
	int num_cell = 0;
	int ele_size = 3;
	int ele_type = 5;
	/** 计算非空气单元数目 **/
	for (int i_tri = 0; i_tri < num_triangle; ++i_tri) {
		int i1 = num_ele - num_triangle + i_tri;
		if (pmeshele[i1].geometry_tag >= tag_air || pmeshele[i1].geometry_tag == tag_remesh) {
			continue;
		}
		num_cell++;
	}

	FILE* fp;
	fp = fopen(fname, "w");
	/*
	 1: points
	 3: line
	 5: Triangular element
	 9: Quadrilateral element
	10: Tetrahedral element
	12: Hexahedral element
	13: Triangular prism element
	14: Pyramid element
	*/
	/** 数据版本声明 **/
	fprintf(fp, "# vtk DataFile Version 2.0\n");
	/** 标题 **/
	fprintf(fp, "vtk title\n");
	/** 文件格式声明 **/
	fprintf(fp, "ASCII\n");
	/** 几何拓扑结构 **/
	fprintf(fp, "DATASET UNSTRUCTURED_GRID\n");
	/** 节点 **/
	fprintf(fp, "\nPOINTS %d double\n", num_pts);
	for (int i = 0; i < num_pts; ++i) {
		fprintf(fp, "%10.8e %10.8e %10.8e\n", pmeshnode[i].x, pmeshnode[i].y, pmeshnode[i].z);
	}
	/** 单元 **/
	fprintf(fp, "\nCELLS %d %d\n", num_cell, (ele_size + 1)*num_cell);
	for (int i_tri = 0; i_tri < num_triangle; ++i_tri) {
		int i1 = num_ele - num_triangle + i_tri;
		if (pmeshele[i1].geometry_tag >= tag_air || pmeshele[i1].geometry_tag == tag_remesh) {
			continue;
		}
		fprintf(fp, "%d", ele_size);
		for (int i = 0; i < ele_size; ++i) {
			fprintf(fp, " %d", pmeshele[i1].n[i]);
		}
		fprintf(fp, "\n");
	}
	/** 单元类型 **/
	fprintf(fp, "\nCELL_TYPES %d\n", num_cell);
	for (int i = 0; i < num_triangle; ++i) {
		int i1 = num_ele - num_triangle + i;
		if (pmeshele[i1].geometry_tag >= tag_air || pmeshele[i1].geometry_tag == tag_remesh) {
			continue;
		}
		fprintf(fp, "%d\n", ele_type);
	}
	/** 数据 **/
	fprintf(fp, "\nPOINT_DATA %d\n", num_pts);
	fprintf(fp, "SCALARS %s double %d\n", "Az/(A/m)", 1);
	fprintf(fp, "LOOKUP_TABLE %s\n", "Atable");
	for (int i = 0; i < num_pts; ++i) {
		if (pmeshnode[i].x > 0) {
			fprintf(fp, "%10.8e\n", A(i) / pmeshnode[i].x);
		} else {
			fprintf(fp, "%10.8e\n", A(i));
		}
	}

	/** 写电磁力 **/
	fprintf(fp, "\nSCALARS %s double %d\n", "F/N", 1);
	fprintf(fp, "LOOKUP_TABLE %s\n", "FTable");
	for (int i = 0; i < num_pts; ++i) {
		double f = 0;
		for (int k = 0; k < 3; k++) {
			f += pmeshnode[i].F[k] * pmeshnode[i].F[k];
		}
		f = sqrt(f);
		fprintf(fp, "%10.8e\n", f);
	}

	/** 写电磁力 **/
	fprintf(fp, "\nVECTORS %s double\n", "FVector");
	for (int i = 0; i < num_pts; ++i) {
		fprintf(fp, "%10.8e %10.8e %10.8e\n", pmeshnode[i].F[0], pmeshnode[i].F[1], pmeshnode[i].F[2]);
	}

	fprintf(fp, "\nCELL_DATA %d\n", num_cell);
	fprintf(fp, "SCALARS %s double %d\n", "B/T", 1);
	fprintf(fp, "LOOKUP_TABLE %s\n", "Btable");
	for (int i_tri = 0; i_tri < num_triangle; ++i_tri) {
		int i1 = num_ele - num_triangle + i_tri;
		if (pmeshele[i1].geometry_tag >= tag_air || pmeshele[i1].geometry_tag == tag_remesh) {
			continue;
		}
		fprintf(fp, "%10.8e\n", pmeshele[i1].B);
	}

	fprintf(fp, "\nVECTORS %s double\n", "Bvector");
	for (int i_tri = 0; i_tri < num_triangle; ++i_tri) {
		int i1 = num_ele - num_triangle + i_tri;
		if (pmeshele[i1].geometry_tag >= tag_air || pmeshele[i1].geometry_tag == tag_remesh) {
			continue;
		}
		fprintf(fp, "%10.8e %10.8e %10.8e\n", pmeshele[i1].Bx, pmeshele[i1].By, 0.0);
	}

	fclose(fp);
}


/*!
 * @brief      计算区域index的电磁力，如果没有指定的话，就按衔铁算。
 * @author     Poofee
 * @date       2020年9月18日
 * @param[out] 
 * @param[in]  int index  
 * @return     void  
**/
double Relay2D::calcMagForce(int index, FORCEMETHOD method) {
    double Fy;
    if(method == FORCEMETHOD::VIRTUALWORK){
        Fy = nodalForceVW(index);
    }
    if(method == FORCEMETHOD::MAXWELLTENSOR){
        Fy = nodalForceMST(index);
    }
    return Fy;
}

/*!
 * @brief      计算二维轴对称电磁力，使用Maxwell's stress tensor
 * 参考
 * 1.Spasov V , Kostov I , Hadzhiev I . Implementation of a
 * novel force computation method in the FEMM software[C]//
 * VII-th International Conference “Engineering, technologies
 * and systems” (TECHSYS 2018), Journal of the Technical
 * University – Sofia, Plovdiv Branch, Bulgaria, "Fundamental
 * Sciences and Applications". 2018.
 * @author     Poofee
 * @date       2020年9月18日
 * @param[out]
 * @param[in]  int index
 * @return     void
**/
double Relay2D::nodalForceMST(int index) {
    if (index < 0) {
        index = tag_xiantie;
    }
    /** 计算运动的节点 **/
    int* moveFlag = (int*)malloc(num_pts * sizeof (int));
    int MOVED = 1;
    for (int i = 0; i < num_pts; ++i) {
        moveFlag[i] = 0;
    }
    for (int i_tri = 0; i_tri < num_triangle; ++i_tri) {
        int k, m, n;
        int i1 = num_ele - num_triangle + i_tri;
        CElement* ele = pmeshele + i1;
        k = ele->n[0]; m = ele->n[1]; n = ele->n[2];
        if (ele->geometry_tag == tag_xiantie) {
            moveFlag[k] = MOVED;
            moveFlag[m] = MOVED;
            moveFlag[n] = MOVED;
        }
    }

    for (int i = 0; i < num_pts; ++i) {
        pmeshnode[i].F[0] = 0.0;
        pmeshnode[i].F[1] = 0.0;
        pmeshnode[i].F[2] = 0.0;
    }
    double Fy = 0, Fx = 0;
    for (int i_tri = 0; i_tri < num_triangle; ++i_tri) {
        int k, m, n;
        int i1 = num_ele - num_triangle + i_tri;
        CElement* ele = pmeshele + i1;
        k = ele->n[0]; m = ele->n[1]; n = ele->n[2];
        CMaterial* mat = getMaterial(ele->geometry_tag);

        /** 寻找空气中的形变单元 **/
        if (ele->geometry_tag == tag_remesh) {
            double rc = ele->rc;
            double mu = mat->getMiu(ele->B);

            double Tensor[2][2];
            Tensor[0][0] = (ele->Bx*ele->Bx-0.5*ele->B*ele->B)/mu;
            Tensor[0][1] = (ele->Bx*ele->By)/mu;
            Tensor[1][0] = (ele->By*ele->Bx)/mu;
            Tensor[1][1] = (ele->By*ele->By-0.5*ele->B*ele->B)/mu;

            for(int j = 0; j < 3; ++j){
                int p = ele->n[j];
                if(moveFlag[p] == MOVED){
                    /** x方向 **/
                    pmeshnode[p].F[0] -= (Tensor[0][0]*ele->P[j]+Tensor[0][1]*ele->Q[j])*PI*rc;
                    /** y方向 **/
                    pmeshnode[p].F[1] -= (Tensor[1][0]*ele->P[j]+Tensor[1][1]*ele->Q[j])*PI*rc;
                }
            }
        }
    }
    for (int i = 0; i < num_pts; ++i) {
        if (moveFlag[i] == MOVED) {
            Fx += pmeshnode[i].F[0];
            Fy += pmeshnode[i].F[1];
        }
    }
    printf("Fx:%10.8e,Fy:%10.8e\n", Fx, Fy);

    free(moveFlag);
    return Fy;
}


/*!
 * @brief      计算二维轴对称电磁力，虚功法
 * 参考
 * 1.Benhama A , Williamson A C , Reece A B J . Virtual work
 * approach to the computation of magnetic force distribution
 * from finite element field solutions[J]. Electric Power
 * Applications, IEE Proceedings -, 2000, 147(6):437-442.
 * 2.戴魏, 余海涛, 胡敏强. 基于虚功法的直线同步电机电磁力计算[J].
 * 中国电机工程学报, 2006, 026(022):110-114.
 * @author     Poofee
 * @date       2020年9月18日
 * @param[out]
 * @param[in]  int index
 * @return     void
**/
double Relay2D::nodalForceVW(int index) {
    if (index < 0) {
        index = tag_xiantie;
    }
    /** 计算运动的节点 **/
    int* moveFlag = (int*)malloc(num_pts * sizeof (int));
    int MOVED = 1;
    for (int i = 0; i < num_pts; ++i) {
        moveFlag[i] = 0;
    }
    for (int i_tri = 0; i_tri < num_triangle; ++i_tri) {
        int k, m, n;
        int i1 = num_ele - num_triangle + i_tri;
        CElement* ele = pmeshele + i1;
        k = ele->n[0]; m = ele->n[1]; n = ele->n[2];
        if (ele->geometry_tag == tag_xiantie) {
            moveFlag[k] = MOVED;
            moveFlag[m] = MOVED;
            moveFlag[n] = MOVED;
        }
    }

    for (int i = 0; i < num_pts; ++i) {
        pmeshnode[i].F[0] = 0.0;
        pmeshnode[i].F[1] = 0.0;
        pmeshnode[i].F[2] = 0.0;
    }
    double Fy = 0, Fx = 0;
    for (int i_tri = 0; i_tri < num_triangle; ++i_tri) {
        int k, m, n;
        int i1 = num_ele - num_triangle + i_tri;
        CElement* ele = pmeshele + i1;
        k = ele->n[0]; m = ele->n[1]; n = ele->n[2];
        CMaterial* mat = getMaterial(ele->geometry_tag);

        /** 寻找空气中的形变单元 **/
        if (ele->geometry_tag == tag_remesh) {
            double rc = ele->rc;
            double area = ele->AREA;
            double mu = mat->getMiu(ele->B);
            double QA = ele->Q[0]*A(k)+ele->Q[1]*A(m)+ele->Q[2]*A(n);
            double PA = ele->P[0]*A(k)+ele->P[1]*A(m)+ele->P[2]*A(n);
            double tmp = PI/4/mu;
            double dvdB2 = mat->getdvdB2(ele->B);
            double tmp1 = PI/4*(QA*QA+PA*PA)/area/rc*dvdB2/4;
            double dQdr[3][3];
            double dPdz[3][3];
            double ac = area*rc;
            double ac2 = ac*ac;
            double ac3 = ac2*ac;

            dQdr[0][0] = 0; dQdr[0][1] = 1; dQdr[0][2] = -1;
            dQdr[1][0] = -1;dQdr[1][1] = 0; dQdr[1][2] = 1;
            dQdr[2][0] = 1; dQdr[2][1] = -1;dQdr[2][2] = 0;

            dPdz[0][0] = 0; dPdz[0][1] = -1;dPdz[0][2] = 1;
            dPdz[1][0] = 1; dPdz[1][1] = 0; dPdz[1][2] = -1;
            dPdz[2][0] = -1;dPdz[2][1] = 1; dPdz[2][2] = 0;

            double dSdr,dSdz;
            double dQAdr,dPAdz;

            for(int j = 0; j < 3; ++j){
                int p = ele->n[j];
                if(moveFlag[p] == MOVED){
                    /** x方向 **/
                    dSdr = 0.5*(dQdr[j][2]*ele->P[1]-dQdr[j][1]*ele->P[2]);
                    dQAdr = (dQdr[j][0]*A(k)+dQdr[j][1]*A(m)+dQdr[j][2]*A(n));
                    pmeshnode[p].F[0] -= tmp*(2*QA/ac*dQAdr);
                    pmeshnode[p].F[0] -= tmp*((PA*PA+QA*QA)*(-1/ac2*(area/3+rc*dSdr)));
                    pmeshnode[p].F[0] -= tmp1*(2*QA/ac2*dQAdr);
                    pmeshnode[p].F[0] -= tmp1*((PA*PA+QA*QA)*(-2/ac3*(area/3+rc*dSdr)));
                    /** y方向 **/
                    dSdz = 0.5*(dPdz[j][1]*ele->Q[2]-dPdz[j][2]*ele->Q[1]);
                    dPAdz = (dPdz[j][0]*A(k)+dPdz[j][1]*A(m)+dPdz[j][2]*A(n));
                    pmeshnode[p].F[1] -= tmp*((PA*PA+QA*QA)*(-1)/ac2*rc*dSdz);
                    pmeshnode[p].F[1] -= tmp*(2*PA/ac*dPAdz);
                    pmeshnode[p].F[1] -= tmp1*((PA*PA+QA*QA)*(-2)/ac3*rc*dSdz);
                    pmeshnode[p].F[1] -= tmp1*(2*PA/ac2*dPAdz);
                }
            }
        }
    }
    for (int i = 0; i < num_pts; ++i) {
        if (moveFlag[i] == MOVED) {
            Fx += pmeshnode[i].F[0];
            Fy += pmeshnode[i].F[1];
        }
    }
    printf("Fx:%10.8e,Fy:%10.8e\n", Fx, Fy);

    free(moveFlag);

    return Fy;
}

/*!
 * @brief      
 * @author     Poofee
 * @date       2020年10月24日
 * @param[out] 
 * @return     void  
**/
void Relay2D::StaticInitial() {
	printf("Please finish it in subclass: void Relay2D::StaticInitial()\n");
}

/*!
 * @brief      
 * @author     Poofee
 * @date       2020年10月24日
 * @param[out] 
 * @return     void  
**/
void Relay2D::TransientInitial() {
	printf("Please finish it in subclass: void Relay2D::TransientInitial()\n");
}

/*!
 * @brief      进行静磁场分析，一般在子类中实现。
 * @author     Poofee
 * @date       2020年10月24日
 * @param[out] 
 * @return     void  
**/
void Relay2D::StaticAnalysis() {
	printf("void Relay2D::StaticAnalysis()\n");
	StaticInitial();
	openGeo();
    MagnetoStatic2DAxisNR();
	/** 输出信息 **/
	outputResults();
}

/*!
 * @brief      进行瞬态分析，一般在子类中实现。
 * @author     Poofee
 * @date       2020年10月24日
 * @param[out] 
 * @return     void  
**/
void Relay2D::TransientAnalysis() {
	printf("void Relay2D::TransientAnalysis()\n");
	TransientInitial();
	openGeo();
    MagnetoDynamics2DAxisNR();
	/** 输出信息 **/
	outputResults();
}
